/*-
 * Copyright (c) 2013 Shang-Rong Cai. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 * $ MTBLAS: stblas_3_gmm.c    Rev. 0.0.20    2013-11-15T12:08+08:00 $
 */

#include <stdlib.h>
//#ifdef TRACE_INFO
#include <stdio.h>
//#endif

#include <bscerr.h>

#include <exec_isa.h>
#include <cor_ctrl.h>
#include <stblas.h>

/* Locally Common Function For Setting Operands of Matrix-Matrix Multiplication */
void stgmm_opchunks_mxmxmulacc(int flag, step_t *step_one, step_t *step_fix, step_t *step_lda, step_t *step_ldb, step_t *step_ldc, addr_t ma, addr_t mb, addr_t mc, size_t c_sz, size_t kz, opchunk_t *pk0a, opchunk_t *pk0b, opchunk_t *pk0c, opchunk_t *pk1a, opchunk_t *pk1b, opchunk_t *pk1c)
{
    pk0a->oc = (addr_t) op4_exec_op;
    pk0a->sz = c_sz;

    pk0a->ra.addr = pk1a;
    pk0a->ra.step = *step_fix;
    pk0a->rb.addr = pk1a;
    pk0a->rb.step = *step_fix;
    pk0a->rc.addr = pk1b;
    pk0a->rc.step = *step_fix;
    pk0a->rd.addr = pk1c;
    pk0a->rd.step = *step_fix;

    pk0b->sz = kz;

    pk0b->ra.addr = ma;
    pk0b->rb.addr = mb;
    pk0b->rb.step = *step_fix;
    pk0b->rc.addr = mc;
    pk0b->rd.addr = NULL;
    pk0b->rd.step = *step_fix;

    pk0c->ra.addr = step_fix;
    pk0c->ra.step = *step_fix;
    pk0c->rb.step = *step_fix;
    pk0c->rc.step = *step_fix;
    pk0c->rd.addr = step_fix;
    pk0c->rd.step = *step_fix;

    pk1a->sz = kz;

    pk1a->ra.addr = ma;
    pk1a->rb.addr = mb;
    pk1a->rc.addr = mc;
    pk1a->rc.step = *step_fix;
    pk1a->rd.addr = NULL;
    pk1a->rd.step = *step_fix;

    pk1b->sz = kz;

    pk1b->ra.addr = ma;
    pk1b->ra.step = *step_fix;
    pk1b->rb.addr = mb;
    pk1b->rc.addr = mc;
    pk1b->rd.addr = NULL;
    pk1b->rd.step = *step_fix;

    pk1c->ra.step = *step_fix;
    pk1c->rb.step = *step_fix;
    pk1c->rc.addr = step_fix;
    pk1c->rc.step = *step_fix;
    pk1c->rd.addr = step_fix;
    pk1c->rd.step = *step_fix;

    if (flag == GEMM_CATEGORY_0)
    {
        pk0b->ra.step = *step_lda;
        pk0b->rc.step = *step_ldc;

        pk0c->rb.addr = step_one;
        pk0c->rc.addr = step_one;

        pk1a->ra.step = *step_one;
        pk1a->rb.step = *step_ldb;

        pk1b->rb.step = *step_one;
        pk1b->rc.step = *step_one;

        pk1c->ra.addr = step_one;
        pk1c->rb.addr = step_ldb;
    }
    else if (flag == GEMM_CATEGORY_1)
    {
        pk0b->ra.step = *step_lda;
        pk0b->rc.step = *step_ldc;

        pk0c->rb.addr = step_ldb;
        pk0c->rc.addr = step_one;

        pk1a->ra.step = *step_one;
        pk1a->rb.step = *step_one;

        pk1b->rb.step = *step_ldb;
        pk1b->rc.step = *step_ldc;

        pk1c->ra.addr = step_one;
        pk1c->rb.addr = step_one;
    }
    else if (flag == GEMM_CATEGORY_2)
    {
        pk0b->ra.step = *step_one;
        pk0b->rc.step = *step_ldc;

        pk0c->rb.addr = step_one;
        pk0c->rc.addr = step_one;

        pk1a->ra.step = *step_lda;
        pk1a->rb.step = *step_ldb;

        pk1b->rb.step = *step_ldb;
        pk1b->rc.step = *step_one;

        pk1c->ra.addr = step_lda;
        pk1c->rb.addr = step_ldb;
    }
    else if (flag == GEMM_CATEGORY_3)
    {
        pk0b->ra.step = *step_one;
        pk0b->rc.step = *step_ldc;

        pk0c->rb.addr = step_ldb;
        pk0c->rc.addr = step_one;

        pk1a->ra.step = *step_lda;
        pk1a->rb.step = *step_one;

        pk1b->rb.step = *step_ldb;
        pk1b->rc.step = *step_one;

        pk1c->ra.addr = step_lda;
        pk1c->rb.addr = step_one;
    }
    else if (flag == GEMM_CATEGORY_4)
    {
        pk0b->ra.step = *step_one;
        pk0b->rc.step = *step_one;

        pk0c->rb.addr = step_ldb;
        pk0c->rc.addr = step_ldc;

        pk1a->ra.step = *step_lda;
        pk1a->rb.step = *step_one;

        pk1b->rb.step = *step_ldb;
        pk1b->rc.step = *step_ldc;

        pk1c->ra.addr = step_lda;
        pk1c->rb.addr = step_one;
    }
    else if (flag == GEMM_CATEGORY_5)
    {
        pk0b->ra.step = *step_one;
        pk0b->rc.step = *step_one;

        pk0c->rb.addr = step_one;
        pk0c->rc.addr = step_ldc;

        pk1a->ra.step = *step_lda;
        pk1a->rb.step = *step_ldb;

        pk1b->rb.step = *step_one;
        pk1b->rc.step = *step_ldc;

        pk1c->ra.addr = step_lda;
        pk1c->rb.addr = step_ldb;
    }
    else if (flag == GEMM_CATEGORY_6)
    {
        pk0b->ra.step = *step_lda;
        pk0b->rc.step = *step_one;

        pk0c->rb.addr = step_ldb;
        pk0c->rc.addr = step_ldc;

        pk1a->ra.step = *step_one;
        pk1a->rb.step = *step_one;

        pk1b->rb.step = *step_ldb;
        pk1b->rc.step = *step_ldc;

        pk1c->ra.addr = step_one;
        pk1c->rb.addr = step_one;
    }
    else if (flag == GEMM_CATEGORY_7)
    {
        pk0b->ra.step = *step_lda;
        pk0b->rc.step = *step_one;

        pk0c->rb.addr = step_one;
        pk0c->rc.addr = step_ldc;

        pk1a->ra.step = *step_one;
        pk1a->rb.step = *step_ldb;

        pk1b->rb.step = *step_one;
        pk1b->rc.step = *step_ldc;

        pk1c->ra.addr = step_one;
        pk1c->rb.addr = step_ldb;
    }
}


void exp_stgmm_opchunks_mxmxmulacc(int flag, step_t *step_one, step_t *step_fix, step_t *step_lda, step_t *step_ldb, step_t *step_buf, addr_t ma, addr_t mb, addr_t buf, size_t c_sz, size_t kz, opchunk_t *pk0a, opchunk_t *pk0b, opchunk_t *pk0c, opchunk_t *pk1a, opchunk_t *pk1b, opchunk_t *pk1c)
{
    pk0a->oc = (addr_t) op4_exec_op;
    pk0a->sz = c_sz;

    pk0a->ra.addr = pk1a;
    pk0a->ra.step = *step_fix;
    pk0a->rb.addr = pk1a;
    pk0a->rb.step = *step_fix;
    pk0a->rc.addr = pk1b;
    pk0a->rc.step = *step_fix;
    pk0a->rd.addr = pk1c;
    pk0a->rd.step = *step_fix;

    pk0b->sz = kz;

    pk0b->ra.addr = ma;
    pk0b->rb.addr = mb;
    pk0b->rb.step = *step_fix;
    pk0b->rc.addr = buf;
    pk0b->rd.addr = NULL;
    pk0b->rd.step = *step_fix;

    pk0c->ra.addr = step_fix;
    pk0c->ra.step = *step_fix;
    pk0c->rb.step = *step_fix;
    pk0c->rc.step = *step_fix;
    pk0c->rd.addr = step_fix;
    pk0c->rd.step = *step_fix;

    pk1a->sz = kz;

    pk1a->ra.addr = ma;
    pk1a->rb.addr = mb;
    pk1a->rc.addr = buf;
    pk1a->rc.step = *step_fix;
    pk1a->rd.addr = NULL;
    pk1a->rd.step = *step_fix;

    pk1b->sz = kz;

    pk1b->ra.addr = ma;
    pk1b->ra.step = *step_fix;
    pk1b->rb.addr = mb;
    pk1b->rc.addr = buf;
    pk1b->rd.addr = NULL;
    pk1b->rd.step = *step_fix;

    pk1c->ra.step = *step_fix;
    pk1c->rb.step = *step_fix;
    pk1c->rc.addr = step_fix;
    pk1c->rc.step = *step_fix;
    pk1c->rd.addr = step_fix;
    pk1c->rd.step = *step_fix;
}


/* Locally Common Function For Setting Operands of Scaled Matrix-Matrix Addition */
void stgmm_opchunks_smma(int rdr, step_t *step_one, step_t *step_fix, step_t *step_lda, step_t *step_buf, addr_t alpha, addr_t ma, addr_t buf, size_t c_sz, opchunk_t *pk_a, opchunk_t *pk_b, opchunk_t *pk_c)
{
    pk_a->sz = c_sz;

    pk_a->ra.addr = alpha;
    pk_a->ra.step = *step_fix;
    pk_a->rb.addr = ma;
    pk_a->rc.addr = buf;
    pk_a->rd.addr = ma;

    pk_b->sz = c_sz;

    pk_b->ra.addr = alpha;
    pk_b->ra.step = *step_fix;
    pk_b->rb.addr = ma;
    pk_b->rc.addr = buf;
    pk_b->rd.addr = ma;

    pk_c->ra.addr = step_fix;
    pk_c->ra.step = *step_fix;
    pk_c->rb.step = *step_fix;
    pk_c->rc.step = *step_fix;
    pk_c->rd.step = *step_fix;

    if (rdr == MTBLAS_ROW_MAJOR)
    {
        pk_a->rb.step = *step_one;
        pk_a->rc.step = *step_one;
        pk_a->rd.step = *step_one;

        pk_b->rb.step = *step_lda;
        pk_b->rc.step = *step_buf;
        pk_b->rd.step = *step_lda;

        pk_c->rb.addr = step_one;
        pk_c->rc.addr = step_one;
        pk_c->rd.addr = step_one;
    }
    else
    {
        pk_a->rb.step = *step_lda;
        pk_a->rc.step = *step_buf;
        pk_a->rd.step = *step_lda;

        pk_b->rb.step = *step_one;
        pk_b->rc.step = *step_one;
        pk_b->rd.step = *step_one;

        pk_c->rb.addr = step_lda;
        pk_c->rc.addr = step_buf;
        pk_c->rd.addr = step_lda;
    }
}


/* Locally Common Function For Setting Operands of Scaled Matrix-Matrix Addition */
void stgmm_opchunks_mxacc(int rdr, step_t *step_one, step_t *step_fix, step_t *step_lda, step_t *step_ldb, addr_t ma, addr_t mb, size_t c_sz, opchunk_t *pk_a, opchunk_t *pk_b, opchunk_t *pk_c)
{
    pk_a->sz = c_sz;

    pk_a->ra.addr = ma;
    pk_a->rb.addr = mb;
    pk_a->rc.addr = NULL;
    pk_a->rc.step = *step_fix;
    pk_a->rd.addr = NULL;
    pk_a->rd.step = *step_fix;

    pk_b->sz = c_sz;

    pk_b->ra.addr = ma;
    pk_b->rb.addr = mb;
    pk_b->rc.addr = NULL;
    pk_b->rc.step = *step_fix;
    pk_b->rd.addr = NULL;
    pk_b->rd.step = *step_fix;

    pk_c->ra.step = *step_fix;
    pk_c->rb.step = *step_fix;
    pk_c->rc.addr = step_fix;
    pk_c->rc.step = *step_fix;
    pk_c->rd.addr = step_fix;
    pk_c->rd.step = *step_fix;

    if (rdr == MTBLAS_ROW_MAJOR)
    {
        pk_a->ra.step = *step_one;
        pk_a->rb.step = *step_one;

        pk_b->ra.step = *step_lda;
        pk_b->rb.step = *step_ldb;

        pk_c->ra.addr = step_one;
        pk_c->rb.addr = step_one;
    }
    else
    {
        pk_a->ra.step = *step_lda;
        pk_a->rb.step = *step_ldb;

        pk_b->ra.step = *step_one;
        pk_b->rb.step = *step_one;

        pk_c->ra.addr = step_lda;
        pk_c->rb.addr = step_ldb;
    }
}


/* Locally Common Function For Setting Operands of Matrix-Matrix Copy */
void stgmm_opchunks_mxmxcopy(int rdr, step_t *step_one, step_t *step_fix, step_t *step_buf, addr_t val, addr_t buf, size_t c_sz, opchunk_t *pk_a, opchunk_t *pk_b, opchunk_t *pk_c)
{
    pk_a->sz = c_sz;
    pk_a->ra.addr = val;
    pk_a->ra.step = *step_fix;
    pk_a->rb.addr = buf;
    pk_a->rc.addr = NULL;
    pk_a->rc.step = *step_fix;
    pk_a->rd.addr = NULL;
    pk_a->rd.step = *step_fix;

    pk_b->sz = c_sz;
    pk_b->ra.addr = val;
    pk_b->ra.step = *step_fix;
    pk_b->rb.addr = buf;
    pk_b->rc.addr = NULL;
    pk_b->rc.step = *step_fix;
    pk_b->rd.addr = NULL;
    pk_b->rd.step = *step_fix;

    pk_c->ra.addr = step_fix;
    pk_c->ra.step = *step_fix;
    pk_c->rb.step = *step_fix;
    pk_c->rc.addr = step_fix;
    pk_c->rc.step = *step_fix;
    pk_c->rd.addr = step_fix;
    pk_c->rd.step = *step_fix;

    if (rdr == MTBLAS_ROW_MAJOR)
    {
        pk_a->rb.step = *step_one;
        pk_b->rb.step = *step_buf;
        pk_c->rb.addr = step_one;
    }
    else
    {
        pk_a->rb.step = *step_buf;
        pk_b->rb.step = *step_one;
        pk_c->rb.addr = step_buf;
    }
}


/* Locally Common Function For Setting Operands of Scaled Matrix */
void stgmm_opchunks_mxscale(int rdr, step_t *step_one, step_t *step_fix, step_t *step_lda, addr_t alpha, addr_t ma, size_t c_sz, opchunk_t *pk_a, opchunk_t *pk_b, opchunk_t *pk_c)
{
#ifdef TRACE_INFO
    printf("[L] stgmm_opchunks_mxscale()\n");
#endif

    pk_a->sz = c_sz;
    pk_a->ra.addr = alpha;
    pk_a->ra.step = *step_fix;
    pk_a->rb.addr = ma;
    pk_a->rc.addr = ma;
    pk_a->rd.addr = NULL;
    pk_a->rd.step = *step_fix;

    pk_b->sz = c_sz;
    pk_b->ra.addr = alpha;
    pk_b->ra.step = *step_fix;
    pk_b->rb.addr = ma;
    pk_b->rc.addr = ma;
    pk_b->rd.addr = NULL;
    pk_b->rd.step = *step_fix;

    pk_c->ra.addr = step_fix;
    pk_c->ra.step = *step_fix;
    pk_c->rb.step = *step_fix;
    pk_c->rc.step = *step_fix;
    pk_c->rd.addr = step_fix;
    pk_c->rd.step = *step_fix;

    if (rdr == MTBLAS_ROW_MAJOR)
    {
        pk_a->rb.step = *step_one;
        pk_a->rc.step = *step_one;
        pk_b->rb.step = *step_lda;
        pk_b->rc.step = *step_lda;
        pk_c->rb.addr = step_one;
        pk_c->rc.addr = step_one;
    }
    else
    {
        pk_a->rb.step = *step_lda;
        pk_a->rc.step = *step_lda;
        pk_b->rb.step = *step_one;
        pk_b->rc.step = *step_one;
        pk_c->rb.addr = step_lda;
        pk_c->rc.addr = step_lda;
    }

#ifdef TRACE_INFO
    printf("    stgmm_opchunks_mxscale() PASSED\n");
#endif
}


int stblas_2e_mxselfscalf(const int rdr, const size_t rz, const size_t cz, const float alpha, const float *ma, const size_t lda)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(float);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(float) * lda;

    const float buf_init = 0;

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_2e_mxselfscalf()\n");
#endif

    if (alpha == 0)
    {
        pk_a.oc = (addr_t) op2_copy_ff;
        pk_b.oc = (addr_t) op2_copy_ff;

        stgmm_opchunks_mxmxcopy(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (addr_t) &buf_init, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

        ra.addr = (addr_t) &pk_a;
        ra.step = step_fix;
        rb.addr = (addr_t) &pk_a;
        rb.step = step_fix;
        rc.addr = (addr_t) &pk_b;
        rc.step = step_fix;
        rd.addr = (addr_t) &pk_c;
        rd.step = step_fix;

        execute((void*) op4_exec_op, rz, &ra, &rb, &rc, &rd);
    }
    else
    {
        if (alpha != 1.0)
        {
            pk_a.oc = (addr_t) op3_mul_ff;
            pk_b.oc = (addr_t) op3_mul_ff;

            stgmm_opchunks_mxscale(rdr, (step_t *) &step_one, (step_t *) &step_fix, (step_t *) &step_lda, (addr_t) &alpha, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_2e_mxselfscalf() PASSED\n");
#endif

    return ierr;
}


int stblas_3e_mxmxmulaccf(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const float *ma, const size_t lda, const float *mb, const size_t ldb, float *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(float);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(float) * lda;
    const step_t step_ldb = sizeof(float) * ldb;
    const step_t step_ldc = sizeof(float) * ldc;

    int flag;  /* Type of matrix order and transpose */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

    opchunk_t pk1a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk1b;  /* Modifier of address */
    opchunk_t pk1c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_3e_mxmxmulaccf()\n");
#endif

    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    pk0b.oc = (addr_t) op3_fmac_ff;
    pk1a.oc = (addr_t) op3_fmac_ff;
    pk1b.oc = (addr_t) op3_fmac_ff;

    stgmm_opchunks_mxmxmulacc(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_ldb, (step_t*) &step_ldc, (addr_t) ma, (addr_t) mb, (addr_t) mc, cz, kz, &pk0a, &pk0b, &pk0c, &pk1a, &pk1b, &pk1c);

    ra.addr = (addr_t) &pk0a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk1b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk0b;
    rc.step = step_fix;
    rd.addr = (addr_t) &pk0c;
    rd.step = step_fix;

    execute(op4_exec_op, rz, &ra, &rb, &rc, &rd);

#ifdef TRACE_INFO
    printf("    stblas_3e_mxmxmulaccf() PASSED\n");
#endif

    return ierr;
}


int stblas_gemmf(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const float alpha, const float *ma, const size_t lda, const float *mb, const size_t ldb, const float beta, float *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(float);
    const step_t step_fix = 0;
    const step_t step_ldc = sizeof(float) * ldc;

    size_t ld_buf;
    size_t step_buf;  /* For buffer matrix leading order */

    int flag;  /* Type of matrix order and transpose */

    float *buf;  /* Buffer vector for storing ma * mb */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_gemmf()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_GEMMF;
    }

    /* Assign buf = mc or Allocate Memory For Buffer: buf[] */
    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        buf = (float*) calloc(rz * cz, sizeof(float));

        if (rdr == MTBLAS_ROW_MAJOR)
        {
            ld_buf = cz;
            step_buf = sizeof(float) * cz;
        }
        else
        {
            ld_buf = rz;
            step_buf = sizeof(float) * rz;
        }

        /* buf[] = 0 */
        ierr = stblas_2e_mxselfscalf(rdr, rz, cz, 0, buf, ld_buf);
    }
    else
    {
        buf = mc;
        ld_buf = ldc;
        step_buf = step_ldc;
    }

    /* mc[] = beta * mc[] */
    ierr = stblas_2e_mxselfscalf(rdr, rz, cz, beta, mc, ldc);

    /* buf[] = ma[] * mb[], buf = alpha * buf[], mc[] = buf[] + mc[] */
    if (alpha != 0)
    {
        /* buf[] = ma[] * mb[] */
        ierr = stblas_3e_mxmxmulaccf(rdr, trna, trnb, rz, cz, kz, ma, lda, mb, ldb, buf, ld_buf);

        /* buf[] = alpha * buf[] */
        ierr = stblas_2e_mxselfscalf(rdr, rz, cz, alpha, buf, ld_buf);

        /* mc[] = buf[] + mc[] */
        if (alpha != 1.0 && beta != 0)
        {
            pk_a.oc = (addr_t) op2_accum_ff;
            pk_b.oc = (addr_t) op2_accum_ff;

            stgmm_opchunks_mxacc(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_buf, (step_t*) &step_ldc, (addr_t) buf, (addr_t) mc, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        free(buf);
    }

    EXIT_STBLAS_GEMMF:
#ifdef TRACE_INFO
    printf("    stblas_gemmf() PASSED\n");
#endif

    return ierr;
}


int stblas_2e_mxselfscal(const int rdr, const size_t rz, const size_t cz, const double alpha, const double *ma, const size_t lda)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(double);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(double) * lda;

    const double buf_init = 0;

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_2e_mxselfscal()\n");
#endif

    if (alpha == 0)
    {
        pk_a.oc = (addr_t) op2_copy_df;
        pk_b.oc = (addr_t) op2_copy_df;

        stgmm_opchunks_mxmxcopy(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (addr_t) &buf_init, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

        ra.addr = (addr_t) &pk_a;
        ra.step = step_fix;
        rb.addr = (addr_t) &pk_a;
        rb.step = step_fix;
        rc.addr = (addr_t) &pk_b;
        rc.step = step_fix;
        rd.addr = (addr_t) &pk_c;
        rd.step = step_fix;

        execute((void*) op4_exec_op, rz, &ra, &rb, &rc, &rd);
    }
    else
    {
        if (alpha != 1.0)
        {
            pk_a.oc = (addr_t) op3_mul_df;
            pk_b.oc = (addr_t) op3_mul_df;

            stgmm_opchunks_mxscale(rdr, (step_t *) &step_one, (step_t *) &step_fix, (step_t *) &step_lda, (addr_t) &alpha, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_2e_mxselfscal() PASSED\n");
#endif

    return ierr;
}



int stblas_3e_mxmxmulacc(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double *ma, const size_t lda, const double *mb, const size_t ldb, double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(double);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(double) * lda;
    const step_t step_ldb = sizeof(double) * ldb;
    const step_t step_ldc = sizeof(double) * ldc;

    int flag;  /* Type of matrix order and transpose */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

    opchunk_t pk1a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk1b;  /* Modifier of address */
    opchunk_t pk1c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_3e_mxmxmulacc()\n");
#endif

    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    pk0b.oc = (addr_t) op3_fmac_df;
    pk1a.oc = (addr_t) op3_fmac_df;
    pk1b.oc = (addr_t) op3_fmac_df;

    stgmm_opchunks_mxmxmulacc(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_ldb, (step_t*) &step_ldc, (addr_t) ma, (addr_t) mb, (addr_t) mc, cz, kz, &pk0a, &pk0b, &pk0c, &pk1a, &pk1b, &pk1c);

    ra.addr = (addr_t) &pk0a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk1b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk0b;
    rc.step = step_fix;
    rd.addr = (addr_t) &pk0c;
    rd.step = step_fix;

    execute(op4_exec_op, rz, &ra, &rb, &rc, &rd);

#ifdef TRACE_INFO
    printf("    stblas_3e_mxmxmulacc() PASSED\n");
#endif

    return ierr;
}


int stblas_gemm(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double alpha, const double *ma, const size_t lda, const double *mb, const size_t ldb, const double beta, double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(double);
    const step_t step_fix = 0;
    const step_t step_ldc = sizeof(double) * ldc;

    size_t ld_buf;
    size_t step_buf;  /* For buffer matrix leading order */

    int flag;  /* Type of matrix order and transpose */

    double *buf;  /* Buffer vector for storing ma * mb */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_gemm()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_GEMM;
    }

    /* Allocate Memory For Buffer: buf[] */
    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        buf = (double*) calloc(rz * cz, sizeof(double));

        if (rdr == MTBLAS_ROW_MAJOR)
        {
            ld_buf = cz;
            step_buf = sizeof(double) * cz;
        }
        else
        {
            ld_buf = rz;
            step_buf = sizeof(double) * rz;
        }

        /* buf[] = 0 */
        ierr = stblas_2e_mxselfscal(rdr, rz, cz, 0, buf, ld_buf);
    }
    else
    {
        buf = mc;
        ld_buf = ldc;
        step_buf = step_ldc;
    }

    /* mc[] = beta * mc[] */
    ierr = stblas_2e_mxselfscal(rdr, rz, cz, beta, mc, ldc);

    /* buf[] = ma[] * mb[], buf = alpha * buf[], mc[] = buf[] + mc[] */
    if (alpha != 0)
    {
        /* buf[] = ma[] * mb[] */
        ierr = stblas_3e_mxmxmulacc(rdr, trna, trnb, rz, cz, kz, ma, lda, mb, ldb, buf, ld_buf);

        /* buf[] = alpha * buf[] */
        ierr = stblas_2e_mxselfscal(rdr, rz, cz, alpha, buf, ld_buf);

        /* mc[] = buf[] + mc[] */
        if (alpha != 1.0 && beta != 0)
        {
            pk_a.oc = (addr_t) op2_accum_df;
            pk_b.oc = (addr_t) op2_accum_df;

            stgmm_opchunks_mxacc(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_buf, (step_t*) &step_ldc, (addr_t) buf, (addr_t) mc, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        free(buf);
    }

    EXIT_STBLAS_GEMM:
#ifdef TRACE_INFO
    printf("    stblas_gemm() PASSED\n");
#endif

    return ierr;
}


int stblas_2e_mxselfscall(const int rdr, const size_t rz, const size_t cz, const long double alpha, const long double *ma, const size_t lda)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(long double);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(long double) * lda;

    const long double buf_init = 0;

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_2e_mxselfscall()\n");
#endif

    if (alpha == 0)
    {
        pk_a.oc = (addr_t) op2_copy_lf;
        pk_b.oc = (addr_t) op2_copy_lf;

        stgmm_opchunks_mxmxcopy(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (addr_t) &buf_init, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

        ra.addr = (addr_t) &pk_a;
        ra.step = step_fix;
        rb.addr = (addr_t) &pk_a;
        rb.step = step_fix;
        rc.addr = (addr_t) &pk_b;
        rc.step = step_fix;
        rd.addr = (addr_t) &pk_c;
        rd.step = step_fix;

        execute((void*) op4_exec_op, rz, &ra, &rb, &rc, &rd);
    }
    else
    {
        if (alpha != 1.0)
        {
            pk_a.oc = (addr_t) op3_mul_lf;
            pk_b.oc = (addr_t) op3_mul_lf;

            stgmm_opchunks_mxscale(rdr, (step_t *) &step_one, (step_t *) &step_fix, (step_t *) &step_lda, (addr_t) &alpha, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_2e_mxselfscall() PASSED\n");
#endif

    return ierr;
}


int stblas_3e_mxmxmulaccl(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const long double *ma, const size_t lda, const long double *mb, const size_t ldb, long double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(long double);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(long double) * lda;
    const step_t step_ldb = sizeof(long double) * ldb;
    const step_t step_ldc = sizeof(long double) * ldc;

    int flag;  /* Type of matrix order and transpose */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

    opchunk_t pk1a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk1b;  /* Modifier of address */
    opchunk_t pk1c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_3e_mxmxmulaccl()\n");
#endif

    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    pk0b.oc = (addr_t) op3_fmac_lf;
    pk1a.oc = (addr_t) op3_fmac_lf;
    pk1b.oc = (addr_t) op3_fmac_lf;

    stgmm_opchunks_mxmxmulacc(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_ldb, (step_t*) &step_ldc, (addr_t) ma, (addr_t) mb, (addr_t) mc, cz, kz, &pk0a, &pk0b, &pk0c, &pk1a, &pk1b, &pk1c);

    ra.addr = (addr_t) &pk0a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk1b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk0b;
    rc.step = step_fix;
    rd.addr = (addr_t) &pk0c;
    rd.step = step_fix;

    execute(op4_exec_op, rz, &ra, &rb, &rc, &rd);

#ifdef TRACE_INFO
    printf("    stblas_3e_mxmxmulaccl() PASSED\n");
#endif

    return ierr;
}


int stblas_gemml(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const long double alpha, const long double *ma, const size_t lda, const long double *mb, const size_t ldb, const long double beta, long double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(long double);
    const step_t step_fix = 0;
    const step_t step_ldc = sizeof(long double) * ldc;

    size_t ld_buf;
    size_t step_buf;  /* For buffer matrix leading order */

    int flag;  /* Type of matrix order and transpose */

    long double *buf;  /* Buffer vector for storing ma * mb */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_gemml()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_GEMML;
    }

    /* Assign buf = mc or Allocate Memory For Buffer: buf[] */
    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        buf = (long double*) calloc(rz * cz, sizeof(long double));

        if (rdr == MTBLAS_ROW_MAJOR)
        {
            ld_buf = cz;
            step_buf = sizeof(long double) * cz;
        }
        else
        {
            ld_buf = rz;
            step_buf = sizeof(long double) * rz;
        }

        /* buf[] = 0 */
        ierr = stblas_2e_mxselfscall(rdr, rz, cz, 0, buf, ld_buf);
    }
    else
    {
        buf = mc;
        ld_buf = ldc;
        step_buf = step_ldc;
    }

    /* mc[] = beta * mc[] */
    ierr = stblas_2e_mxselfscall(rdr, rz, cz, beta, mc, ldc);

    /* buf[] = ma[] * mb[], buf = alpha * buf[], mc[] = buf[] + mc[] */
    if (alpha != 0)
    {
        /* buf[] = ma[] * mb[] */
        ierr = stblas_3e_mxmxmulaccl(rdr, trna, trnb, rz, cz, kz, ma, lda, mb, ldb, buf, ld_buf);

        /* buf[] = alpha * buf[] */
        ierr = stblas_2e_mxselfscall(rdr, rz, cz, alpha, buf, ld_buf);

        /* mc[] = buf[] + mc[] */
        if (alpha != 1.0 && beta != 0)
        {
            pk_a.oc = (addr_t) op2_accum_lf;
            pk_b.oc = (addr_t) op2_accum_lf;

            stgmm_opchunks_mxacc(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_buf, (step_t*) &step_ldc, (addr_t) buf, (addr_t) mc, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        free(buf);
    }

    EXIT_STBLAS_GEMML:
#ifdef TRACE_INFO
    printf("    stblas_gemml() PASSED\n");
#endif

    return ierr;
}


int stblas_2e_cmxselfscalf(const int rdr, const size_t rz, const size_t cz, const float complex alpha, const float complex *ma, const size_t lda)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(float complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(float complex) * lda;

    const float complex buf_init = 0;

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_2e_cmxselfscalf()\n");
#endif

    if (alpha == 0)
    {
        pk_a.oc = (addr_t) op2_copy_fc;
        pk_b.oc = (addr_t) op2_copy_fc;

        stgmm_opchunks_mxmxcopy(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (addr_t) &buf_init, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

        ra.addr = (addr_t) &pk_a;
        ra.step = step_fix;
        rb.addr = (addr_t) &pk_a;
        rb.step = step_fix;
        rc.addr = (addr_t) &pk_b;
        rc.step = step_fix;
        rd.addr = (addr_t) &pk_c;
        rd.step = step_fix;

        execute((void*) op4_exec_op, rz, &ra, &rb, &rc, &rd);
    }
    else
    {
        if (alpha != 1.0)
        {
            pk_a.oc = (addr_t) op3_mul_fc;
            pk_b.oc = (addr_t) op3_mul_fc;

            stgmm_opchunks_mxscale(rdr, (step_t *) &step_one, (step_t *) &step_fix, (step_t *) &step_lda, (addr_t) &alpha, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_2e_cmxselfscalf()\n");
#endif

    return ierr;
}


int stblas_3e_cmxmxmulaccf(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const float complex *ma, const size_t lda, const float complex *mb, const size_t ldb, float complex *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(float complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(float complex) * lda;
    const step_t step_ldb = sizeof(float complex) * ldb;
    const step_t step_ldc = sizeof(float complex) * ldc;

    int flag;  /* Type of matrix order and transpose */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

    opchunk_t pk1a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk1b;  /* Modifier of address */
    opchunk_t pk1c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_3e_cmxmxmulaccf()\n");
#endif

    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (trna != MTBLAS_TRANS_CONJ && trnb != MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_fmac_fc;
        pk1a.oc = (addr_t) op3_fmac_fc;
        pk1b.oc = (addr_t) op3_fmac_fc;
    }
    else if (trna == MTBLAS_TRANS_CONJ && trnb != MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_conjfmac_fc;
        pk1a.oc = (addr_t) op3_conjfmac_fc;
        pk1b.oc = (addr_t) op3_conjfmac_fc;
    }
    else if (trna != MTBLAS_TRANS_CONJ && trnb == MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_conjbfmac_fc;
        pk1a.oc = (addr_t) op3_conjbfmac_fc;
        pk1b.oc = (addr_t) op3_conjbfmac_fc;
    }
    else
    {
        pk0b.oc = (addr_t) op3_conjabfmac_fc;
        pk1a.oc = (addr_t) op3_conjabfmac_fc;
        pk1b.oc = (addr_t) op3_conjabfmac_fc;
    }

    stgmm_opchunks_mxmxmulacc(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_ldb, (step_t*) &step_ldc, (addr_t) ma, (addr_t) mb, (addr_t) mc, cz, kz, &pk0a, &pk0b, &pk0c, &pk1a, &pk1b, &pk1c);

    ra.addr = (addr_t) &pk0a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk1b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk0b;
    rc.step = step_fix;
    rd.addr = (addr_t) &pk0c;
    rd.step = step_fix;

    execute(op4_exec_op, rz, &ra, &rb, &rc, &rd);

#ifdef TRACE_INFO
    printf("    stblas_3e_cmxmxmulaccf() PASSED\n");
#endif

    return ierr;
}


int stblas_cgemmf(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const float complex alpha, const float complex *ma, const size_t lda, const float complex *mb, const size_t ldb, const float complex beta, float complex *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(float complex);
    const step_t step_fix = 0;
    const step_t step_ldc = sizeof(float complex) * ldc;

    size_t ld_buf;
    size_t step_buf;  /* For buffer matrix leading order */

    int flag;  /* Type of matrix order and transpose */

    float complex *buf;  /* Buffer vector for storing ma * mb */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_cgemmf()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_CGEMMF;
    }

    /* Assign buf = mc or Allocate Memory For Buffer: buf[] */
    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        buf = (float complex*) calloc(rz * cz, sizeof(float complex));

        if (rdr == MTBLAS_ROW_MAJOR)
        {
            ld_buf = cz;
            step_buf = sizeof(float complex) * cz;
        }
        else
        {
            ld_buf = rz;
            step_buf = sizeof(float complex) * rz;
        }

        /* buf[] = 0 */
        ierr = stblas_2e_cmxselfscalf(rdr, rz, cz, 0, buf, ld_buf);
    }
    else
    {
        buf = mc;
        ld_buf = ldc;
        step_buf = step_ldc;
    }

    /* mc[] = beta * mc[] */
    ierr = stblas_2e_cmxselfscalf(rdr, rz, cz, beta, mc, ldc);

    /* buf[] = ma[] * mb[], buf = alpha * buf[], mc[] = buf[] + mc[] */
    if (alpha != 0)
    {
        /* buf[] = ma[] * mb[] */
        ierr = stblas_3e_cmxmxmulaccf(rdr, trna, trnb, rz, cz, kz, ma, lda, mb, ldb, buf, ld_buf);

        /* buf[] = alpha * buf[] */
        ierr = stblas_2e_cmxselfscalf(rdr, rz, cz, alpha, buf, ld_buf);

        /* mc[] = buf[] + mc[] */
        if (alpha != 1.0 && beta != 0)
        {
            pk0a.oc = (addr_t) op2_accum_fc;
            pk0b.oc = (addr_t) op2_accum_fc;

            stgmm_opchunks_mxacc(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_buf, (step_t*) &step_ldc, (addr_t) buf, (addr_t) mc, cz, &pk0a, &pk0b, &pk0c);

            ra.addr = (addr_t) &pk0a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk0b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk0c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        free(buf);
    }

    EXIT_STBLAS_CGEMMF:
#ifdef TRACE_INFO
    printf("    stblas_cgemmf() PASSED\n");
#endif

    return ierr;
}


int stblas_2e_cmxselfscal(const int rdr, const size_t rz, const size_t cz, const double complex alpha, const double complex *ma, const size_t lda)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(double complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(double complex) * lda;

    const double complex buf_init = 0;

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_2e_cmxselfscal()\n");
#endif

    if (alpha == 0)
    {
        pk_a.oc = (addr_t) op2_copy_dc;
        pk_b.oc = (addr_t) op2_copy_dc;

        stgmm_opchunks_mxmxcopy(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (addr_t) &buf_init, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

        ra.addr = (addr_t) &pk_a;
        ra.step = step_fix;
        rb.addr = (addr_t) &pk_a;
        rb.step = step_fix;
        rc.addr = (addr_t) &pk_b;
        rc.step = step_fix;
        rd.addr = (addr_t) &pk_c;
        rd.step = step_fix;

        execute((void*) op4_exec_op, rz, &ra, &rb, &rc, &rd);
    }
    else
    {
        if (alpha != 1.0)
        {
            pk_a.oc = (addr_t) op3_mul_dc;
            pk_b.oc = (addr_t) op3_mul_dc;

            stgmm_opchunks_mxscale(rdr, (step_t *) &step_one, (step_t *) &step_fix, (step_t *) &step_lda, (addr_t) &alpha, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_2e_cmxselfscal()\n");
#endif

    return ierr;
}


int stblas_3e_cmxmxmulacc(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double complex *ma, const size_t lda, const double complex *mb, const size_t ldb, double complex *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(double complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(double complex) * lda;
    const step_t step_ldb = sizeof(double complex) * ldb;
    const step_t step_ldc = sizeof(double complex) * ldc;

    int flag;  /* Type of matrix order and transpose */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

    opchunk_t pk1a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk1b;  /* Modifier of address */
    opchunk_t pk1c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_3e_cmxmxmulacc()\n");
#endif

    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (trna != MTBLAS_TRANS_CONJ && trnb != MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_fmac_dc;
        pk1a.oc = (addr_t) op3_fmac_dc;
        pk1b.oc = (addr_t) op3_fmac_dc;
    }
    else if (trna == MTBLAS_TRANS_CONJ && trnb != MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_conjfmac_dc;
        pk1a.oc = (addr_t) op3_conjfmac_dc;
        pk1b.oc = (addr_t) op3_conjfmac_dc;
    }
    else if (trna != MTBLAS_TRANS_CONJ && trnb == MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_conjbfmac_dc;
        pk1a.oc = (addr_t) op3_conjbfmac_dc;
        pk1b.oc = (addr_t) op3_conjbfmac_dc;
    }
    else
    {
        pk0b.oc = (addr_t) op3_conjabfmac_dc;
        pk1a.oc = (addr_t) op3_conjabfmac_dc;
        pk1b.oc = (addr_t) op3_conjabfmac_dc;
    }


    stgmm_opchunks_mxmxmulacc(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_ldb, (step_t*) &step_ldc, (addr_t) ma, (addr_t) mb, (addr_t) mc, cz, kz, &pk0a, &pk0b, &pk0c, &pk1a, &pk1b, &pk1c);

    ra.addr = (addr_t) &pk0a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk1b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk0b;
    rc.step = step_fix;
    rd.addr = (addr_t) &pk0c;
    rd.step = step_fix;

    execute(op4_exec_op, rz, &ra, &rb, &rc, &rd);

#ifdef TRACE_INFO
    printf("    stblas_3e_cmxmxmulacc() PASSED\n");
#endif

    return ierr;
}


int stblas_cgemm(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double complex alpha, const double complex *ma, const size_t lda, const double complex *mb, const size_t ldb, const double complex beta, double complex *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(double complex);
    const step_t step_fix = 0;
    const step_t step_ldc = sizeof(double complex) * ldc;

    size_t ld_buf;
    size_t step_buf;  /* For buffer matrix leading order */

    int flag;  /* Type of matrix order and transpose */

    double complex *buf;  /* Buffer vector for storing ma * mb */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_cgemm()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_CGEMM;
    }

    /* Assign buf = mc or Allocate Memory For Buffer: buf[] */
    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        buf = (double complex*) calloc(rz * cz, sizeof(double complex));

        if (rdr == MTBLAS_ROW_MAJOR)
        {
            ld_buf = cz;
            step_buf = sizeof(double complex) * cz;
        }
        else
        {
            ld_buf = rz;
            step_buf = sizeof(double complex) * rz;
        }

        /* buf[] = 0 */
        ierr = stblas_2e_cmxselfscal(rdr, rz, cz, 0, buf, ld_buf);
    }
    else
    {
        buf = mc;
        ld_buf = ldc;
        step_buf = step_ldc;
    }

    /* mc[] = beta * mc[] */
    ierr = stblas_2e_cmxselfscal(rdr, rz, cz, beta, mc, ldc);

    /* buf[] = ma[] * mb[], buf = alpha * buf[], mc[] = buf[] + mc[] */
    if (alpha != 0)
    {
        /* buf[] = ma[] * mb[] */
        ierr = stblas_3e_cmxmxmulacc(rdr, trna, trnb, rz, cz, kz, ma, lda, mb, ldb, buf, ld_buf);

        /* buf[] = alpha * buf[] */
        ierr = stblas_2e_cmxselfscal(rdr, rz, cz, alpha, buf, ld_buf);

        /* mc[] = buf[] + mc[] */
        if (alpha != 1.0 && beta != 0)
        {
            pk0a.oc = (addr_t) op2_accum_dc;
            pk0b.oc = (addr_t) op2_accum_dc;

            stgmm_opchunks_mxacc(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_buf, (step_t*) &step_ldc, (addr_t) buf, (addr_t) mc, cz, &pk0a, &pk0b, &pk0c);

            ra.addr = (addr_t) &pk0a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk0b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk0c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        free(buf);
    }

    EXIT_STBLAS_CGEMM:
#ifdef TRACE_INFO
    printf("    stblas_cgemm() PASSED\n");
#endif

    return ierr;
}


int stblas_2e_cmxselfscall(const int rdr, const size_t rz, const size_t cz, const long double complex alpha, const long double complex *ma, const size_t lda)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(long double complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(long double complex) * lda;

    const long double complex buf_init = 0;

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_2e_cmxselfscall()\n");
#endif

    if (alpha == 0)
    {
        pk_a.oc = (addr_t) op2_copy_lc;
        pk_b.oc = (addr_t) op2_copy_lc;

        stgmm_opchunks_mxmxcopy(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (addr_t) &buf_init, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

        ra.addr = (addr_t) &pk_a;
        ra.step = step_fix;
        rb.addr = (addr_t) &pk_a;
        rb.step = step_fix;
        rc.addr = (addr_t) &pk_b;
        rc.step = step_fix;
        rd.addr = (addr_t) &pk_c;
        rd.step = step_fix;

        execute((void*) op4_exec_op, rz, &ra, &rb, &rc, &rd);
    }
    else
    {
        if (alpha != 1.0)
        {
            pk_a.oc = (addr_t) op3_mul_lc;
            pk_b.oc = (addr_t) op3_mul_lc;

            stgmm_opchunks_mxscale(rdr, (step_t *) &step_one, (step_t *) &step_fix, (step_t *) &step_lda, (addr_t) &alpha, (addr_t) ma, cz, &pk_a, &pk_b, &pk_c);

            ra.addr = (addr_t) &pk_a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk_b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk_c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_2e_cmxselfscall()\n");
#endif

    return ierr;
}


int stblas_3e_cmxmxmulaccl(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const long double complex *ma, const size_t lda, const long double complex *mb, const size_t ldb, double long complex *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(long double complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(long double complex) * lda;
    const step_t step_ldb = sizeof(long double complex) * ldb;
    const step_t step_ldc = sizeof(long double complex) * ldc;

    int flag;  /* Type of matrix order and transpose */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */
    operand_t rd;  /* Usual operand or point to modifier of step */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

    opchunk_t pk1a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk1b;  /* Modifier of address */
    opchunk_t pk1c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_3e_cmxmxmulaccl()\n");
#endif

    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (trna != MTBLAS_TRANS_CONJ && trnb != MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_fmac_lc;
        pk1a.oc = (addr_t) op3_fmac_lc;
        pk1b.oc = (addr_t) op3_fmac_lc;
    }
    else if (trna == MTBLAS_TRANS_CONJ && trnb != MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_conjfmac_lc;
        pk1a.oc = (addr_t) op3_conjfmac_lc;
        pk1b.oc = (addr_t) op3_conjfmac_lc;
    }
    else if (trna != MTBLAS_TRANS_CONJ && trnb == MTBLAS_TRANS_CONJ)
    {
        pk0b.oc = (addr_t) op3_conjbfmac_lc;
        pk1a.oc = (addr_t) op3_conjbfmac_lc;
        pk1b.oc = (addr_t) op3_conjbfmac_lc;
    }
    else
    {
        pk0b.oc = (addr_t) op3_conjabfmac_lc;
        pk1a.oc = (addr_t) op3_conjabfmac_lc;
        pk1b.oc = (addr_t) op3_conjabfmac_lc;
    }

    stgmm_opchunks_mxmxmulacc(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_ldb, (step_t*) &step_ldc, (addr_t) ma, (addr_t) mb, (addr_t) mc, cz, kz, &pk0a, &pk0b, &pk0c, &pk1a, &pk1b, &pk1c);

    ra.addr = (addr_t) &pk0a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk1b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk0b;
    rc.step = step_fix;
    rd.addr = (addr_t) &pk0c;
    rd.step = step_fix;

    execute(op4_exec_op, rz, &ra, &rb, &rc, &rd);

#ifdef TRACE_INFO
    printf("    stblas_3e_cmxmxmulaccl() PASSED\n");
#endif

    return ierr;
}


int stblas_cgemml(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const long double complex alpha, const long double complex *ma, const size_t lda, const long double complex *mb, const size_t ldb, const long double complex beta, long double complex *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(long double complex);
    const step_t step_fix = 0;
    const step_t step_ldc = sizeof(long double complex) * ldc;

    size_t ld_buf;
    size_t step_buf;  /* For buffer matrix leading order */

    int flag;  /* Type of matrix order and transpose */

    long double complex *buf;  /* Buffer vector for storing ma * mb */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to opchunk to target of modifying */
    operand_t rc;  /* Usual operand or point to modifier of address */

    opchunk_t pk0a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk0b;  /* Modifier of address */
    opchunk_t pk0c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_cgemml()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemm_flag(rdr, trna, trnb, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_CGEMML;
    }

    /* Assign buf = mc or Allocate Memory For Buffer: buf[] */
    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        buf = (long double complex*) calloc(rz * cz, sizeof(long double complex));

        if (rdr == MTBLAS_ROW_MAJOR)
        {
            ld_buf = cz;
            step_buf = sizeof(long double complex) * cz;
        }
        else
        {
            ld_buf = rz;
            step_buf = sizeof(long double complex) * rz;
        }

        /* buf[] = 0 */
        ierr = stblas_2e_cmxselfscall(rdr, rz, cz, 0, buf, ld_buf);
    }
    else
    {
        buf = mc;
        ld_buf = ldc;
        step_buf = step_ldc;
    }

    /* mc[] = beta * mc[] */
    ierr = stblas_2e_cmxselfscall(rdr, rz, cz, beta, mc, ldc);

    /* buf[] = ma[] * mb[], buf = alpha * buf[], mc[] = buf[] + mc[] */
    if (alpha != 0)
    {
        /* buf[] = ma[] * mb[] */
        ierr = stblas_3e_cmxmxmulaccl(rdr, trna, trnb, rz, cz, kz, ma, lda, mb, ldb, buf, ld_buf);

        /* buf[] = alpha * buf[] */
        ierr = stblas_2e_cmxselfscall(rdr, rz, cz, alpha, buf, ld_buf);

        /* mc[] = buf[] + mc[] */
        if (alpha != 1.0 && beta != 0)
        {
            pk0a.oc = (addr_t) op2_accum_lc;
            pk0b.oc = (addr_t) op2_accum_lc;

            stgmm_opchunks_mxacc(rdr, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_buf, (step_t*) &step_ldc, (addr_t) buf, (addr_t) mc, cz, &pk0a, &pk0b, &pk0c);

            ra.addr = (addr_t) &pk0a;
            ra.step = step_fix;
            rb.addr = (addr_t) &pk0b;
            rb.step = step_fix;
            rc.addr = (addr_t) &pk0c;
            rc.step = step_fix;

            execute((void*) op4_exec_op, rz, &ra, &ra, &rb, &rc);
        }
    }

    if (alpha != 0 && alpha != 1.0 && beta != 0)
    {
        free(buf);
    }

    EXIT_STBLAS_CGEMML:
#ifdef TRACE_INFO
    printf("    stblas_cgemml() PASSED\n");
#endif

    return ierr;
}
